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We study the quantum phase diagram of a three dimensional non-interacting Dirac semimetal 
in the presence of either quenched axial or scalar potential disorder, by calculating the average 
and the typical density of states as well as the inverse participation ratio using numerically exact 
methods. We show that as a function of the disorder strength a half-filled (i.e. undoped) Dirac 
semimetal displays three distinct ground states, namely an incompressible semimetal, a compressible 
diffusive metal, and a localized Anderson insulator, in stark contrast to a conventional dirty metal 
that only supports the latter two phases. We establish the existence of two distinct quantum 
critical points, which respectively govern the semimetal-metal and the metal-insulator quantum 
phase transitions and also reveal their underlying multifractal nature. Away from half-filling the 
(doped) system behaves as a diffusive metal that can undergo Anderson localization only, which is 
shown by determining the mobility edge and the phase diagram in terms of energy and disorder. 
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The development of Lorentz invariant relativistic quan¬ 
tum mechanics naturally led to the discovery of the cel¬ 
ebrated Dirac equation [lj to describe the electron and 
positron. In high energy physics there are no massless 
Dirac fermions, as both leptons and quarks are massive 
particles and it is only at energy scales higher than the 
corresponding mass gap, that they can be approximated 
as massless. In contrast, many condensed matter sys¬ 
tems support a gapless phase, where two-fold Kramers 
degenerate conduction and valence bands touch linearly 
at isolated points in the Brillouin zone. The low energy 
excitations around these diabolic points (known as Dirac 
points, where the conical conduction and valence bands 
touch each other) can be described by a massless Dirac 
equation in the infrared limit, such a phase is known as 
a Dirac semimetal (DSM) if the system is undoped with 
the Fermi level precisely at the Dirac point. 

Recently, an intense experimental investigation into 
narrow gap semiconductors [1] has led to the discov¬ 
ery of three dimensional DSMs in various materials 
such as Cd 3 As 2 (§- 0 , Na 3 Bi [ 3 , 0 , Bii x Sb x 00 , 
BiTl(Si_ 5 Se 5 )9 10, El, (Bii_ x In x ) 2 Se 3 [l| Q, and 
Pbi_ x Sn x Te |2|, IldL 1151 . These experiments have now 
raised the exciting prospect of studying the physics of 
three dimensional massless Dirac fermions in solid state 
systems. Due to the ubiquitous presence of disorder in 
all solid state systems, a precise understanding of the 
phase diagram of dirty DSMs is therefore a problem of 
deep fundamental importance, which can only be studied 
in condensed matter systems since the corresponding rel¬ 
ativistic Dirac problem in particle physics does not, by 
definition, have any disorder. While effects of disorder 
in conventional metals are fairly well understood HfHUi 
this is not the case for three dimensional DSMs where the 
applicable quantum phase diagram and (even) the ques¬ 


tion of how many allowed phases may exist as a function 
of disorder are still wide open. 


Due to the finite density of states (DOS) at the Fermi 
energy, a conventional metal (CM) (or a Fermi liquid) is a 
compressible state of matter, and an infinitesimally weak 
disorder acts as a relevant perturbation for the ballistic 
Fermi liquid fixed point, giving rise to a finite lifetime for 
the ballistic single particle excitations (20j. Therefore, a 
dirty, noninteracting CM in three dimensions only sup¬ 
ports two phases, namely a diffusive metal and a localized 
Anderson insulator (AI) [2lj. The Anderson localiza¬ 
tion is a continuous quantum phase transition (QPT), for 
which the spatial dimension d = 2 serves as the lower crit¬ 
ical dimension [5lj |. This QPT is reflected in the spatial 
variation of the wave function and can only be captured 
by quantities that probe its extended or localized nature, 
such as the inverse participation ratio (IPR) [3, 231, and 


the typical DOS (TDOS) [H, 3, 25|. However, a self av¬ 
eraging quantity such as the average DOS (ADOS) itself 
remains unaffected by the localization transition. 


In contrast to CMs, a half-filled (i.e. undoped) three 
dimensional DSM is an incompressible fermionic quan¬ 
tum critical system, which leads to a quadratically van¬ 
ishing DOS at the Fermi level. Consequently, any weak 
disorder acts as an irrelevant perturbation for a three 
dimensional DSM. Therefore, the three dimensional bal¬ 
listic DSM phase should remain stable up to a critical 
strength of disorder W c . At this threshold, the DSM un¬ 
dergoes a continuous QPT into a compressible diffusive 
metal (CDM) phase H3|, and the ADOS at zero en¬ 
ergy acts as an order parameter for describing this tran¬ 
sition, while distinguishing it from a conventional An¬ 
derson localization. Based on a nonlinear sigma model 
analysis, the CDM phase has been argued to display an 
Anderson localization at a higher strength of disorder Wt , 
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which belongs to the unitary (A) universality class [2t 
The first analysis of the DSM-CDM transition in Ref. 2£ 
was a mean-field one, and only recently the non-Gaussian 
nature of this quantum critical point (QCP) has been elu¬ 
cidated 271. Subsequently, the DSM-CDM phase tran¬ 
sition has been addressed u sing various analytic [28-31 [ 
and numerical methods 32|-34 1. However, the effects of 
large disorder and the Anderson localization transition 
have not been studied through any numerically exact ap¬ 
proach. Moreover, there have been recent claims that 
rare-region effects convert the DSM into a CDM for an 
infinitesimally weak disorder strength Ebl so that there 
is in fact no DSM-CDM transition at all! Therefore, 
it is still an important open question to figure out how 
many distinct quantum phases indeed exist in a dirty 3D 
Dirac system for finite disorder (one or two or three with 
DSM/CDM/AI being all the possibilities) through a non- 
perturbative calculation, in addition to establishing the 
appropriate universality of any applicable disorder-driven 
QPT in the system. In the current paper we answer this 
question through extensive exact numerical work supple¬ 
mented by theoretical arguments. 




FIG. 1: Phase diagram of three dimensional Dirac semimetals. 
(Left) pa (0) and pt(0) for L = 60, clearly displaying the exis¬ 
tence of three distinct phases DSM, CDM, and AI. The TDOS 
tracks the ADOS inside the DSM, through the semimetal- 
metal QCP, and then pt( 0) goes to zero at Wi, signaling the 
Anderson localization transition. (Right) M e d ge (E) as a func¬ 
tion of energy and disorder strength for L = 60. The white 
and the blue regions are respectively metallic and localized, 
and the blank region is gapped (outside the bandwidth) which 
clearly show the existence of the energy dependent mobility 
edge. The value of M e d ge (E) is shown in the key. 


We study the phase diagram of a dirty DSM in three 
dimensions using numerically exact methods. We calcu¬ 
late the ADOS and the TDOS on sufficiently large lat¬ 
tice sizes reaching up to 60 3 sites, using the kernel poly¬ 
nomial method (KPM) [3(1]. As the TDOS tracks the 
ADOS in any metallic phase and also serves as an order 
parameter for Anderson localization [T^i EH > it can nat¬ 
urally capture both possible (DSM-CDM-AI) QPTs (see 
Fig. CD. We firmly establish that a disordered DSM in 
three dimensions possesses two distinct QPTs as a func¬ 
tion of disorder strength, as shown in Figs. [0 and [3] The 
TDOS at zero energy is only finite in the CDM region 
W c < W < Wi, in contrast to the ADOS at zero energy 


which is finite for any W > W c . This provides unambigu¬ 
ous evidence for the the existence of three phases and the 
different nature of the two distinct disorder-tuned QPTs. 
The nature of the QCP between the DSM and the CDM 
phases has been studied through an extensive numerical 
computation of the ADOS in Refs. 32, 34 ]. As the TDOS 
and the ADOS track each other in a metal, our calcula¬ 
tion provides further numerical evidence for the stability 
of the DSM phase in the presence of weak disorder. We 
supplement the DOS analysis through exact calculations 
of the wave function for small system sizes in order to 
estimate the critical exponents of the Anderson localiza¬ 
tion transition. Here, our goal is not to establish the 
precise values for the critical exponents of the Anderson 
localization transition. Rather, we want to show that 
the estimated critical exponents for the CDM-AI QCP 
are comparable to the ones known for the conventional 
orthogonal (AI) Wigner-Dyson universality class EM, 
and are strikingly different from the ones obtained for the 
DSM-CDM QCP [HEl. Since our model of a DSM pre¬ 
serves time reversal symmetry, we expect the Anderson 
localization transition to be described by the orthogonal 
(AI) class 17, lH, 371 for short range disorder. 


We consider the following Hamiltonian on a cubic lat¬ 
tice with periodic boundary conditions, 


H = \ cvV’r+e* + H.c) + v (r)$A w i/> r 

Z r,£ 

(1) 

where ipj = (Cr^f, c r ,+ 4 , c ri _ 4 ) is a four compo¬ 

nent Dirac spinor composed of an electron at site r, with 
a parity (±) and spin (f / \), jx = x,y,z, ep is a unit 
vector that points to each nearest neighbor, and the a^ 
are 4x4 anticommuting Dirac matrices. In the absence 
of disorder, the tight binding model gives rise to eight 
Dirac cones at the T, M, R and X points of the cubic 
Brillouin zone. The disorder potential at each site V (r) 
is a randomly distributed variable between [—W/2, W/2] 
and Aw is a matrix that specifies the type of disorder. 
We primarily focus on the random axial chemical poten¬ 
tial Aw = 75 = iaia 2 (X 3 1 but also discuss the effects of a 
random scalar potential Aw = ^ 4 x 4 (where Rxi is a four 
by four identity matrix). We consider these two types 
of disorder because we want to demonstrate universal¬ 
ity in the sense that both disorders give rise to identical 
phase diagrams. This is also implied by the existence of 
two continuous U(l) symmetries under ip —> e ie ip and 
ip —> e l ^ lh ip. In addition, axial disorder is pertinent for 
addressing the phase diagram of dirty gapless supercon¬ 
ductors. 

TDOS and ADOS: To track the nature of each phase 
of the model defined in Eq. CD. we first calculate the 
ADOS, the local DOS (LDOS), and the TDOS using the 


























3 


KPM (see Ref. [36| for details), which are defined as 


fashion as 


Pa(E) = 

\ i=1 a=l / 

Pia(E) = El( fc ’ P\ii a )\ 2 K E ~ E k p), 
k,/3 

/ ^ Ns 4 


( 2 ) 

( 3 ) 


p t ( E) = exp ( — EE( log (i;)) ) . (4) 


i =1 ot —1 


The size of the system considered is K = L 3 , where | i, a) 
denotes an eigenstate at site i and orbital a (one of the 
four Dirac orbitals) with an energy Ei a , and the (...) 
denotes a disorder average. As the ADOS is self aver¬ 
aging, we only include the disorder average to smoothen 
the data. For the calculations of the ADOS presented 
here with a system size V = 60 3 we have only used 12 
disorder averages. In contrast, the TDOS is not a self 
averaging quantity and it requires a significant number 
of disorder averages to obtain a convergent result. After 
the disorder averaging, translation symmetry is restored 
and consequently we also sum over a finite number of lat¬ 
tice sites N s -C V to improve the statistics. We provide 
the parameters used in the study in the Supplemental 
Material 37|]. Central to the KPM, we first expand the 
LDOS or the ADOS in terms of the orthogonal Cheby- 
shev polynomials up to an order N c . For the ADOS we 
have considered N c = 1028, whereas the computation of 
the TDOS is much more demanding [HI and for most of 
the calculations we have used N c = 8192, unless stated 
otherwise. 



FIG. 2: Critical properties of the two quantum critical points. 
(Left) Power law dependence of p a ( 0) and pt(0) as a function 
of disorder strength in the vicinity of the DSM-CDM QCP 
for L = 60. (Right) The Anderson localization QCP from the 
TDOS as a function of the Chebyshev expansion order N c for 
L = 30 respectively in three and (Inset) two dimensions. We 
find the localization transition in three dimensions is converg¬ 
ing to Wf d /t = 8.9 ± 0.3 in three dimensions whereas in two 
dimensions W 2d /t —¥ 0 as N c —> oo by extrapolating pt(0) to 
zero (dashed lines). 


The DSM is characterized by an ADOS p a {E) oc \E\ 2 , 
that vanishes at zero energy. As shown in Fig. [TJ both 
p a (0) and p t (0) vanish for weak disorder strengths, and 
concomitantly become finite after passing through the 
DSM-CDM QCP at W c /t « 2.55 [34j], in a power law 


Pa (E = 0, W > W c ) 
Pt {E = 0 1 W>W c ) 


(w - w c y 
(w - w c y 


( 5 ) 

( 6 ) 


where we find x = 1.4 ± 0.2 (in agreement with Refs. [32, 
34l for much larger system sizes) and xt = 2.0 ± 0.3 as 
shown in Fig. |T| Physically, the average moments of the 
LDOS ( Pi{E) n ) captures the multifractal nature of the 
disordered wavefunction EE EE [Hi ■ Therefore the dif¬ 
ference Xt — x = 0.60 ± 0.36 reveals the underlying multi¬ 
fractal nature of the DSM-CDM QCP. We have checked 
that the value of W c determined from p a ( 0) is within 
numerical accuracy unaffected by increasing N c . 

Inside the CDM phase p t { 0) tracks p a (0) {pt ~ p a ) up 
to W/t « 4.5. As p a is a self averaging quantity, we 
can conclude that the DSM-CDM transition completely 
independent of Anderson localization. For larger values 
of disorder ( W/t > 4.5) p a remains finite and p t goes to 
zero continuously at the localization transition as 


Pt(E = 0, W < Wi) ~ {Wi - W) p . 


( 7 ) 


We determine the location of the transition by studying 
the effect of increasing the Chebyshev expansion order on 
p t , as Wi has been shown to be sensitive to N c in Ref. l39l. 
By extrapolating p t (E = 0) to zero, we find Wi/t = 
8.9 ± 0.3 with an order parameter critical exponent /3 = 
1.5 ± 0.2, as shown in Fig. [2] 




FIG. 3: TDOS with L = 60 as a function of energy for vari¬ 
ous values of disorder. (Left) Passing through the DSM-CDM 
QCP, (Right) passing through the Anderson localization tran¬ 
sition. The arrows denote increasing disorder strengths from 
W/t = 2.0 to 4.5 (Left) and from W/t = 5.0 to 13.0 (Right). 
We see no signs of rare region effects [35j in the numerics. 


In order to further understand the stability of the CDM 
against localization effects, we compare the TDOS for a 
DSM in two and three spatial dimensions in Fig. [2] The 
two dimensional model of a DSM is obtained by setting 
the hopping along the 0 direction to zero in Eq. ©, 
which gives rise to four Dirac cones at the high sym¬ 
metry points of the Brillouin zone for a square lattice. 
In two dimensions, both axial and scalar potential dis¬ 
orders are marginally relevant perturbations, and an in¬ 
finitesimal amount of disorder leads to an AI 37l 4oL 411. 
This is in excellent agreement with our numerical results, 
namely W 2d /t —> 0 as N c —> oo. In contrast, in three di¬ 
mensions we find the DSM phase to be stable up to W c 
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(Fig. [0, where p a ( 0) and p t ( 0) simultaneously develop 
finite values, and the CDM phase undergoes an Ander¬ 
son localization at a much larger disorder strength W). 
Therefore, our results establish the existence of three dis¬ 
tinct ground states for a dirty DSM in three dimensions 
(DSM/CDM/AI with increasing disorder), in agreement 
with the field theoretic analysis. 

By considering the energy dependence of pt as shown 
in Fig. [3] we find pt(E ) ~ \E\ 2 inside the DSM phase 
(i.e. in agreement with p a (E ) S3)- In contrast, p t ( 0) 
approaches a constant inside the CDM phase. Within 
this framework, we can capture the location of the mo¬ 
bility edge separating the CDM and AI phases by con¬ 
sidering the ratio M e d ge (E) = p t (E)/p a (E). As shown 
in Fig. [T| starting from the CDM phase, the system de¬ 
velops a finite, energy dependent mobility edge and the 
entire band localizes for a larger disorder strength going 
to an AI phase, qualitatively similar to the situation for a 
three dimensional CM. Due to the underlying Dirac band 
structure, the shape of the mobility edge is different from 
a conventional tight binding model (e.g. see Ref. HU). 
This phase diagram also demonstrates that a DSM away 
from half-filling behaves as a CM. 

Wave function: In order to understand the critical 
properties of the localization transition and compare our 
results with the well established universality classes for 
Anderson localization, we now study the properties of the 
wave function. The qualitative nature of the wave func¬ 
tion also displays all of the physics we have discussed so 
far [37]. Here, we focus on the average participation ratio 
defined as 


P — 

- 1 avg — 


212 


Ki,J^MrQl 2 ] 

Ei.a IVu(u )| 4 


( 8 ) 


for wave functions i/s a ( r i) and only focus on the center 
of the band E = 0 (otherwise Pf ih, will be a function 
of energy). It is well known |l8| that the IPR scales 
as P-Jg ~ 1/V in a metallic phase and P~J- ~ const, 
in a localized phase [37]. As these calculations re¬ 
quire an exact diagonalization of the Hamiltonian, we 
are restricted to much smaller lattice sizes than consid¬ 
ered for the DOS. Here we focus on linear system sizes 
L = 4,6,8,10,12 with 10,000 disorder realizations for 
L = 4—10 and 1,000 for L = 12. Due to the limited 
system sizes we first expand the IPR in the vicinity of 
the localization transition in terms of an unknown scal¬ 
ing function / and L dependent corrections, [II Elfin 
P-\ « L~ d2 (/((1 - Wi/W)L x ' v ) + A 0 /L y ) . We have 
introduced the fractal dimension c^, the localization 
length exponent is, is a finite size correction, and y the 
leading order irrelevant variable’s critical exponent fiiSj ]. 
To determine the location of the transition we calculate 
the the best recursive fit to Yipr = P a AL d2 — A 0 L~ V , us¬ 
ing finite size scaling techniques 43j, [46§, see Fig. [IJ As a 
result, we find a localization transition at Wi/t = 8.8±0.3 




FIG. 4: Finite size scaling of the IPR. (Left) Determining the 
location of the Anderson localization transition by optimizing 
the crossing location of Yipr by third order polynomial fits to 
the IPR, this yields Wi/t = 8.8±0.3. (Right) Data collapse of 
Yipr in the vicinity of Wi as a function of (1 — Wi/W)L 1 ^ 1 ', 
yields a localization length exponent v = 1.5 ± 0.2. (Inset) 
Optimizing the data collapse using the local linearity function 
S(v) (as described in Ref. [48j), the value of v is determined 
by the minimum of S(is). 


in excellent agreement with the typical DOS results. This 
also yields c ?2 = 1.2±0.3, y = 2.3±0.5, and Aq = 1.8±0.6. 

Now that we have determined the critical point us¬ 
ing two different methods we perform scaling data col¬ 
lapse to extract the localization length exponent, which 
yields is = 1.5 ± 0.2. Within the numerical accuracy, 
v , y, and d 2 are consistent with the known exponents 
for the orthogonal (AI) class (is = 1.57, y = 2.8, and 
c ?2 = 1.23) |42|, |43[ [47[. In previous studies of Ander¬ 
son localization, it was found that typical DOS exponent 
satisfies /3 = is(aa — d ) Em. We also find this rela¬ 
tion to hold with op = 4.0 ± 0.3, which is consistent with 
ao = 4.043 fi3 for class AI . 


In conclusion, we have performed a detailed numeri¬ 
cally exact study of the quantum phase diagram and An¬ 
derson localization in three dimensional dirty DSMs. By 
studying the TDOS on sufficiently large lattices we have 
established the existence of three phases for a half-filled 
system DSM, CDM and AI, separated by two different 
QCPs as a function of disorder. We have supplemented 
this with a precise finite size scaling analysis of the IPR 
and obtained the critical exponents of the localization 
transition. We also establish clear signatures for rnulti- 
fractility at the QPTs. In addition, we find the localiza¬ 
tion properties for both axial and scalar potential disor¬ 
ders to be identical (to within numerically accuracy) fiH , 
and therefore, we conclude that the QPTs driven by ei¬ 
ther disorder belong to the same universality class. As 
the model we have studied for axial or potential disor¬ 
der separates into two blocks of Weyl Hamiltonians |37j, 
our results are equally valid for three dimensional time 
reversal symmetric Weyl systems. 
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Supplemental Material 

In the supplemental material we first provide field theoretic arguments to justify our claims regarding the universality 
class of the localization transition. We then describe in detail the parameters used in the numerical study, present the 
typical DOS for axial and potential disorder, the system size dependence of p t { 0), the finite size scaling of the IPR 
for potential disorder, present a comparison between the typical and average density of states in each phase studied 
in this work, and show both the inverse participation ratio for each system size studied and the wave function in each 
phase as a function of the distance to its maximal absolute value. 


FIELD THEORY ARGUMENTS 

In momentum space, the tight-binding Dirac Hamiltonian (Eq. (1) in the main text) in the absence of disorder 
becomes 


3 

H = t sin k^a^ (SI) 

M=1 


with the dispersion relation 


E k = ±t y sin 2 k x + sin 2 k y + sin 2 k z . (S2) 

The conduction and the valence bands touch at following eight high symmetry points (i) (0,0,0), (ii) (n, 7r, 0), (iii) 
(7r, 0,7r), (iv) (0,7r, 7r), (v) (7r, 0, 0), (vi) (0,7r,0), (vii) (0,0,7r), and (viii) (7r,7r,7r) of the cubic Brillouin zone, as a 
consequence of the Fermion doubling theorem. The lattice Hamiltonian commutes with 75 = iaia 2 oi 3 , which generates 
a (7(1) chiral symmetry. Due to the absence of any mass term, the Hamiltonian naturally anti-commutes with two 
other Dirac matrices /3 and /J 75 . The linearized Dirac Hamiltonian around these points can be written as 

H D j = Vijk x ai + v 2 ,jkya 2 + v 3t] k z a 3 . (S3) 

The velocities (v\,v 2 ,v 3 ) at the eight Dirac points are respectively given by (v,v,v), (—v,—v,v), (—v,v,—v), 
(v,—v,—v), (—v,v,v), (v,—v,v), (v,v,—v) and (—v,—v,—v), where v = t (in units of H = a = 1 for a lattice 
spacing a). With suitable unitary transformations at different cones we can covert all of the linearized Hamiltonians 
at cones (ii) to (viii) into the same form as that of (i). For example the spinor ipg at (n, 7 t, 7 t) can be transformed 
as /3ipg or ifi'ygipg, to arrive at the desired form. Thus, the linearized theory will have an emergent SU( 8 ) flavor 
symmetry. This is a continuum approximation to the cubic point group symmetry. Similarly, if we set the hopping 
in the z direction to zero, in order to obtain the two dimensional Hamiltonian studied in the main text, we obtain 
four species of four-component Dirac fermions in the vicinity of (0,0), (w, 0), (0, 7 r) and ( 7 r, 7 r). The linearized theory 
then demonstrates an SU(4:) flavor symmetry. Any inter-valley scattering process due to disorder and/or interaction 
usually reduces this emergent symmetry. 

The discussion of disorder effects can be easily understood in the chiral basis, with 07 = ay ® 73, a 2 = cr 2 <8> t 3 , 
a 3 = 03 <B> 73, /3 = <Jo ti and 75 = <8) 73. Here, 07 are Pauli matrices acting on the spin index, and 17 are 
Pauli matrices which act on the chiral indices. In this basis the lattice Hamiltonian in the presence of both axial and 
potential disorders is block diagonal and consists of the following two 2x2 Weyl Hamiltonians 

Hl,W = (^Xpr^Xl.r+eA + H.c) + ^[Kx(r) + ^pWjxLx'l.r, 

r 

H2,W = ~ (^X2,r °>X 2 ,r+e A + H - c ) + + V p( r )]X2,rX2,r, 

r,fi r 

where Xi/ 2 r = ( c r, 1 / 2 ,u c r, 1/2.4.) a ^ w0 component spinor composed of the linear combinations c r y i(T = [c r ,+,a + 

c r -^\/y/2 and c r , 2 ,a = [cr,+,tr — Cr,-,(j]/\/2- Each block by itself represents a model Hamiltonian of a Weyl semimetal 
on a non-centrosymmetric cubic lattice, which has four right handed and four left handed Weyl points [49(. For this 
reason axial and potential disorders show identical behavior, and the model clearly has time reversal symmetry. For 
each two component lattice model (in the presence of the intervalley scattering due to short range disorder), the 


(54) 

(55) 
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relevant symmetry class is the orthogonal Wigner-Dyson class or AI. Consequently, the CDM to AI transition will 
belong to this class. Similar to the case of graphene with inter-valley scattering |40|, |4l|, in the replica formalism the 
appropriate Q matrix for this problem will be chosen from the space sp{ 2 n)^Sp( 2 n) • the P resence of intervalley 
scattering (which is always a relevant term inside the diffusive metal phase) the Q matrices for different Weyl cones 
are locked into a single one. In Ref. [2(1 the Q matrix was chosen to be a unitary matrix, which is actually applicable 
for a time reversal symmetry breaking Weyl semimetal. It is also interesting to recall that the identical behaviors 
of the axial and the scalar potential disorder for the DSM-CDM transition were first noted in the Supplementary 
Material of Ref. [27]- The renormalization group flow with respect to the ballistic DSM fixed point in the presence of 
both types of disorder was found to be 


jA l 

dl 


A j + 2 A 


2 

V 


z - 1 +Ab 

i 


(5 6 ) 

(57) 


where z is the scale dependent dynamic exponent, and and A p denote the dimensionless coupling constants for 
axial and potential disorder respectively. The DSM-CDM transition occurs along A CjJ 4 + A Cj p = 1/2, with z c \ = 3/2 
and t'd = 1. This further supports our symmetry based argument for the identical behaviors of the scalar potential 
and the axial disorders, which has also been verified through extensive numerical calculations j34j- Note that we are 
using the subscripts 1 and 2 to denote the DSM-CDM and CDM-AI quantum phase transitions respectively. 

For completeness we also write down the following nonlinear sigma model for the 4-component Dirac Hamiltonian, 

S[Q] = J d 3 x Tr(VQ) 2 + J d 3 x Tr[HQ], (S8) 

where p a (0) ~ § v ci(.d-z c i) ; g ^he avera g e DOS inside the CDM phase, D is the diffusion constant inside the CDM phase 
at the zero energy and 5 = 2(A a + Ay — 1/2) is the reduced distance from the DSM-CDM critical point. The product 
p a (0)D(0) ~ <Jxx{ 0) ~ 5 Vcl ^~ d ~ 2 \ where a xx is the longitudinal conductivity of the CDM. In the above equation D is 
the external frequency. The dynamic critical exponent for the Anderson localization is z c 2 = d = 3, which is clearly 
distinct from z c 1 =3/2 for the DSM-CDM QCP. The beta function for this model is well known up to five loop order 
and can be found in Refs. IHoMHij 


ADDITIONAL NUMERICAL RESULTS 


TABLE SI: Table of the parameters used for the typical DOS. We denote the number of disorder realizations used as Nd and 
the number of lattice sites that we average the typical DOS over is N s it e . 


L 

N d 

N s ite 

20 

2000 

20 

30 

800 

10 

40 

700 

10 

50 

600 

10 

60 

520 

10 


The overall behavior of the wave function also captures the localization transition, while it is unaffected by the 
DSM-CDM QCP, see Fig. [SU We compute the decay of the wave function from its maximal value ’F(r) = |t/(|r — 
r max|)|/|^( r max)| for the center of the band, where r max is the location in real space where |^(r)| is maximal. We 
find no real qualitative change in the wave function across the DSM-CDM QCP, in complete agreement with the fact 
that this is not a localization transition in the wave function. In contrast, we find a drastic change in the qualitative 
nature of the wave function across the CDM-AI localization transition, where for W > Wi , 'P(r) ~ exp(— |r —r max |/A) 
with a finite localization length A in the AI phase. 




















FIG. SI: (Left) The typical density of states at zero energy as a function of disorder W with L = 20 for axial and potential 
disorder. Remarkably, we find the results are equivalent to within numerical accuracy and therefore conclude the critical 
exponent /3 is the same for these two types of disorder. (Right) Typical density of states for various system sizes, which show 
a very weak dependence on the system size for L > 30. 
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FIG. S2: Finite size scaling of the IPR for potential disorder. (Left) Determining the location of the Anderson localization using 
third order polynomial fits to obtain Wi = 8.4 ±0.3, which agrees with the TDOS within numerical accuracy. We attribute the 
small finite shift from the TDOS result for Wi to finite size effects. From this we obtain the critical exponents d 2 = 1-2 ± 0.3, 
y = 2.0 ± 0.5, Ao = 1.8 ± 0.6. (Right) Data collapse of the IPR, which yields v = 1.6 ± 0.2, in excellent agreement with the 
axial disorder result. (Inset) Quality of the scaling collapse is dictated by minimizing the local linearity function. 
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FIG. S3: Comparison of the typical and average density of states with L = 60 in the Dirac semi-metal phase W/t = 2.0 (a), 
the diffusive metal phase W/t = 3.5 (b) and W/t = 5.0 (c), and the Anderson insulating phase W/t = 12.5 (d). We see the 
typical density of states tracks the average in both the Dirac phase and diffusive metal phase and goes to zero for all energies 
in the Anderson insulating phase. 



W/t 


FIG. S4: Inverse participation ratio versus disorder strength for each system size studied. We find the expected behavior in 
a delocalized phase (pertaining to both the DSM and CDM phases) ~ 1/V and in the localized phase approaching a 
constant. 
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FIG. S5: Decay of the magnitude of the wave function 'F(r) = |-0(|r — rmaKDI/I^Crmax)! for the center of the band as a function 
of the distance to its maximal value |'i/ , (r m ax)| for L = 10 in each phase studied. We find no real qualitative change in the wave 
function across the DSM to CDM QCP (left and center), whereas in the Anderson insulating phase the wave function displays 
a clear decay (right). 

















